Topological persistence and dynamical heterogeneities near jamming 
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We introduce topological methods for quantifying spatially heterogeneous dynamics, and use these 
tools to analyze particle-tracking data for a quasi-two-dimensional granular system of air-fluidized 
beads on approach to jamming. In particular we define two overlap order parameters, which quantify 
the correlation between particle configurations at different times, based on a Voronoi construction 
and the persistence in the resulting cells and nearest neighbors. Temporal fluctuations in the 
decay of the persistent area and bond order parameters define two alternative dynamic four-point 
susceptibilities, Xa(j) and Xb(t~), well-suited for characterizing spatially-heterogeneous dynamics. 
These are analogous to the standard four-point dynamic susceptibility Xi{^ r )> but where the space- 
dependence is fixed uniquely by topology rather than by discretionary choice of cutoff function. 
While these three susceptibilities yield characteristic time scales that are somewhat different, they 
give domain sizes for the dynamical heterogeneities that are in good agreement and that diverge on 
approach to jamming. 

PACS numbers: 45.70.n, 47.55.Lm, 64.70.Pf, 02.40.Pc 



I. INTRODUCTION 

The dynamics in both thermal and driven systems be- 
come increasingly heterogeneous [l|, H, Q on approach 
to jamming, where a control parameter such as temper- 
ature, packing fraction, or driving rate is varied to the 
point that all rearrangements cease [3, [f| . Such spatially- 
heterogeneous dynamics (SHD) may take the form of in- 
termittent strings or swirls of neighboring particles that 
follow one another in correlated fashion from one pack- 
ing configuration to another. Observations of SHD have 
been reported for coarsening @,J3 and sheared [1] foams, 
for dense colloidal suspensions [jj [l(| HH > an d f° r driven 
granular systems [HI, [r|. Quantification of SHD may 
be accomplished by a dynamic four-point susceptibility, 
Xi(J-: T )i defined in terms of density correlations across 
intervals of both space and time [14| . Bounds on XiQ: T ) 
have been inferred from standard measurements of di- 
electric susceptibility for glass-forming liquids and of dy- 
namic light scattering for dense colloidal suspensions [15| . 
Direct measurements of Xi(^ T ) have been reported for 
quasi-two-dimensional granular systems, where particle 
motion was excited by shear [16] and by a fluidizing up- 
flow of air [171 ]. In the former, Xi(^ T ) was measured vs 
space and time at a fixed packing fraction; in the lat- 
ter, Xi{li T ) was measured vs time at fixed I for several 
packing fractions. In both Refs. [HI [l3], the dynamic 
susceptibilities are based on temporal fluctuations of an 
overlap order parameter; however, the former employs 
a Gaussian cutoff function while the latter uses a step 
function. 

In this paper we conduct further measurements and 
analyses of spatially-heterogeneous dynamics (SHD) in 
a system of air-fluidized beads similar to Ref. \u\ . Our 
primary extension of Ref. (l7| concerns the cutoff func- 
tion in the overlap order parameter. First we explore the 
length-dependence of Xi(h T ) by varying the cutoff dis- 
tance I. Second, we remove the arbitrariness in choice 



of cutoff function by using topological measures of over- 
lap. In particular we build upon the concept of persis- 
tence, introduced to help characterize the coarsening of 
spin domains in magnets or of bubbles in foams (l8l . Il9| . 
As a topological measure of overlap for particulate sys- 
tems, we define a normalized, dimcnsionless persistent 
area as the fraction of space that remains inside the same 
Voronoi cell after a given time interval. We also introduce 
a second topological measure of configurational overlap 
in terms of persistent bonds, i.e. the fraction of nearest- 
neighbor pairs that remain neighbors after a given time 
interval. While both persistent area and persistent bond 
order parameters vanish with time, the latter can exhibit 
a much slower decay if a string of particles follow one an- 
other over a long distance. After discussion of such cutoff 
functions, we use the resulting four-point susceptibilities 
both to explore the connection of SHD to local structure 
and to observe the evolution of SHD on approach to jam- 
ming. Finally we consider possible artifacts due to finite 
experiment duration. 



II. METHODS AND BACKGROUND 

The granular system under study here is a quasi- 
two-dimensional monolayer of spheres that roll with- 
out slipping due to a fluidizing upfiow of air, as in 
Refs. Ho, [HI- The spheres are steel with an equal 
number of d s = 0.318 cm and di = 0.397 cm diameter 
sizes. 



This is the same system as in Ref. 17], though 
here most run durations are 120 minutes rather than 
just 20 minutes. The steel beads are confined in a cir- 
cular flat sieve with a diameter of 17.7 cm and a mesh 
size of 150 /im. Energy is injected by a uniform flow 
of air through the sieve at a flux which is high enough 
to randomly drive the balls by turbulence (Re 10 3 ) 
without causing levitation. The flux is increased from 
560 cm/s to 700 cm/s as the area fraction is decreased 



FIG. 1: (Color online) (a-c) Example average velocity vector 
fields, Ar/r, for time intervals r and bead area fractions (j) 
as specified in (d). At short times, r < Tb as in (a), the 
bead motion is ballistic and there are no spatial correlations. 
At very long times, as in (c), the bead motion is diffusive 
and there are no spatial correlations. At intermediate times, 
t c < t < T r , the bead motion is sub-diffusive with beads 
remaining temporarily trapped in a cage of fixed neighbors. 
Toward the end of this caging regime, as in (b), the bead 
motion exhibits spatially heterogeneous dynamics in the form 
of string-like swirls in the average velocity vector field. This 
figure was generated from data described in Ref. [2(J, where 
the bead diameters are 0.873 cm and 0.635 cm, and where 
precise definitions of time-scales are given in terms of the 
logarithmic slope of the mean-squared displacement vs time. 
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FIG. 2: (Color online) Time dependence of (a) the aver- 
age overlap order parameter and (b) the corresponding four- 
point susceptibilities, defined by Eqs. (12I3[) respectively, for 
air-fluidized beads at an area packing fraction of <j> — 0.792. 
Different curves are for different cutoff lengths, I, as labeled. 
The thick green curve is for the standard cutoff threshold 
I — 0.5d s given by half the small-bead diameter. 
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FIG. 3: (Color online) Time dependence of (a) the average 
overlap order parameter and (b) the corresponding four-point 
susceptibilities, for air-fluidized beads at various area frac- 
tions, as labeled. The cutoff threshold is I — 0.5d s . 



from (f> = 0.8 to (j) = 0.6, to ensure that the low density 
systems are fluidized. Three feet above the sieve are six 
incandescent lights aligned in a ring one foot in diam- 
eter. At the center of the ring is a 120 Hz Pulnix 6710 
video camera with a 644 x 484 array of 8BIT square pixels 
and a zoom lens. Light reflects specularly off the tops of 
the beads and is imaged by the camera. Bright spots are 
tracked so that the positions r$ (t) are known for all beads 
i = {1, 2, . . . , N} and for all times t during the entire ex- 
periment. Further details are available in Refs. [20l. [2l|. 

By increasing the fraction <f> of area occupied by 
spheres, at fixed gas flux, this system exhibits a jamming 
transition at point-J, such that the average bead kinetic 
energy vanishes linearly on approach to random close 
packing at about 4> c = 0.83 [2(J. Near this transition, the 
bead dynamics exhibit the usual sequence of crossovers 
from ballistic, to sub-diffusive, to diffusive motion as a 
function of delay time. Spatially-heterogeneous dynam- 
ics occur in the sub-diffusive regime, particularly at delay 
times such that the beads are breaking out of their cages 
and beginning to diffuse. This is illustrated graphically 
in Fig. [TJa-c) , which show example average velocity vec- 
tor fields, Ar/r, for three different delay times r; the 
area fraction and delay time for each image is chosen to 
highlight motion in each of the three regions of phase 
space, as specified in Fig. [TJi. For both very short and 
very long delay times, corresponding to ballistic and dif- 
fusive regimes, the magnitude and direction of the aver- 
age velocities for neighboring beads are completely uncor- 
rected. However, at intermediate delay times, Fig. [TH, 
corresponding to subdiffusive motion, a few string-like 
clusters of neighboring beads move in swirls while other 
beads move relatively little. As time evolves, such swirls 
come and go in different regions of the sample so that the 
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dynamics become completely ergodic. This is best seen 
by a video of the average velocity field [22j . 

Now we begin the central task of this paper, which is to 
quantify the spatially-heterogenous nature of the dynam- 
ics depicted qualitatively in Fig. [T] For this, a standard 
tool is the four-point dynamic susceptibility X4(^ T )i as 
reviewed in Ref. [lij]. This may be approximated from 
temporal fluctuations in the instantaneous self-overlap 
order parameter, defined here as 



dynamical correlation length, £, which diverges at jam- 
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where N is the number of beads and Wi is a cutoff func- 
tion that equals 1 if the displacement of bead i remains 
less than I across the whole time interval t — > t + r and 
that equals otherwise. Whereas here and in Ref. [l7| 
Wi is a step-function cutoff, in Ref. [HI ] it is a Gaussian. 
In either case the average self-overlap order parameter 
and the four-point susceptibility are given by moments 
of Qt(l, t) averaged over all times t as follows: 
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Note that X4(^ T ) is independent of system size because 
the variance scales as 1/N. Also note that X4 grows in 
proportion to the variance, or heterogeneity, in the dis- 
placements of different beads across a given time interval. 

Example results for the space and time dependence of 
Q(1,t) and X4G? r ) are shown in Fig. [2J for one partic- 
ular area fraction <j> — 0.792. In particular, these two 
functions are plotted vs delay time r for several different 
values of I. By construction Q(1,t) decays from one to 
zero as a function of r. For increasing I, the observed 
decay becomes somewhat sharper, but more importantly 
the location of the decay moves to longer delay times, 
since the beads take longer to move a larger distance. 
By contrast the susceptibility X4^> T ) is a peaked func- 
tion of t, which vanishes at both early and late times. 
By comparison of Figs, d^-b it is evident that the peak 
of \4(1,t) is located at the characteristic decay time of 
Q{1,t). Furthermore, the height of the peak depends on 
I, such that it vanishes for large and small I and becomes 
a maximum for I on the order of the bead size. Thus 
in Ref. [13] the overlap threshold was set equal to the 
radius of the small beads, I = 0.5d s , and the dynamic 
susceptibility was considered only as a function of time, 
Xa(t) = X4(0.5d s ,r). 

Results for Q(t) and Xi( T ) vs T are displayed in 
Figs. [3^,-b for a sequence of area fractions. This repro- 
duces Fig. 3c of Ref. [13], with the identical system but 
with a longer run duration. As seen before, the location 
t| of the peak in X4 (r) moves to longer times on approach 
to jamming. Furthermore, the peak height xt increases 
concurrently. Thus the dynamics become not only slower 
but more spatially heterogeneous on approach to jam- 
ming. The size of the spatial heterogeneities represents a 



ming; such behavior was first reported in [17| . and will 
be discussed further in Section IV. Note that for all area 
fractions, the value of the overlap order parameter is 
nearly constant, Q* = 0.48 ± 0.06, when the suscepti- 
bility reaches its peak. 

The dynamical correlation length may be deduced from 
the value of the self-overlap order parameter and the 
dynamic susceptibility at its peak, using the following 
physical picture. For a system of N total beads, we 
wish to know the typical number M of heterogeneities 
and more importantly the typical number n of beads 
in each of these fast-moving domains. Since the con- 
tribution to the self-overlap order parameter is Qq = 
from the nM highly-mobile beads inside domains, and 
is Qi = 1 from the other N — nM less-mobile beads, 
the average self-overlap order parameter at T4 is Q* = 
(QM)) = [Qo[nM) + Q X {N - nM)]/N = 1 - nM/N. 
However the instantaneous value of Qt(T|) varies with 
time t because of counting statistics AM — \J~M in the 
number of heterogeneous domains actually present at any 
instant. Therefore, the peak dynamic susceptibility is 
X\ = N(AQ*) 2 = n 2 M/N. By eliminating M/N in these 
expressions for Q* and x&> t ne typical number of parti- 
cles in each fast-moving dynamic heterogeneity is found 
to be 



714 



1 



(4) 



Since Q* is nearly constant, as observed in Fig. [3l the 
number of beads per domain is directly proportional to 
the peak susceptibility x\- If the domains are com- 
pact, then the dynamical correlation length scales as 
S,/d s ~ n x l d 1 where d s is the bead diameter and d is di- 
mensionality [IBj[l3|- ^ the domains are strings, then the 
dynamical correlation length is £/d s ~ n. If fluctuations 
An in domain size are also included, then the left-hand 
side of Eq. ([JJ is slightly modified to n[l + (An/n) 2 ]. 

To our knowledge, the counting arguments leading to 
Eq. 0$ have not been previously published. While it 
seems common knowledge that the peak susceptibility 
depends upon domain size, we have found no arguments 
or calculations in prior literature to support or quantify 
the relation. 



III. TOPOLOGICAL PERSISTENCE 

The dynamic susceptibility X4(^ r ) discussed in the 
previous section relies upon a somewhat arbitrary choice 
for a cutoff function, w, in Eq. {J). Furthermore, since 
t is the more interesting parameter, it also relies upon 
a particular choice of cutoff length, I, so that spatially- 
heterogeneous dynamics may be studied as a function of 
delay time. In this section we propose two alternative dy- 
namic susceptibilities, intended for the same purpose, in 
which the role of the cutoff function and cutoff length are 
determined uniquely by topological considerations. Both 
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are based upon the Voronoi tessellation constructed from 
the particle positions at each instant of time. 



A. Persistent area 

An order parameter for quantifying the degree of over- 
lap between particle configurations may be constructed 
by considering the fraction of the sample which remains 
inside the same Voronoi cell across a time interval t — » 
t + r. In particular, we define the instantaneous persis- 
tent volume as 



Mt) = V J at ^ T ^ dv 



(5) 



where V is the sample volume, a t (r,r) = 1 if the point 
r remains inside the same Voronoi cell across the whole 
time interval t — ► t + r, and <it(r, r) = if the point r 
becomes enclosed by another Voronoi cell. Since our ex- 
periments are in two-dimensions, we shall refer to Eq. ([5|) 
as the instantaneous persistent area. By construction, it 
decays monotonically from one to zero as the beads move 
a distance comparable to their size. Thus the persistent 
area is similar to the previous self-overlap order parame- 
ter Qt{l, t) with a cutoff length I set by bead size. 

A graphic illustration of the contribution a t (r, r) of 
each point in the sample to the average over space given 
by Eq. ([5]) is shown in Fig. 3] for a sequence of increasing 
delay times r. Pixels are colored black for at iT (r) = 1 
and are colored white otherwise; therefore, the sample is 
all black for r = and it becomes progressively white 
for increasing r. As illustrated in Fig. [4j white regions 
first appear at short r near the boundary of the initial 
Voronoi cells. At larger r the white regions thicken as the 
persistent black regions at the cell cores shrink. Eventu- 
ally each black spot vanishes. The graphics here can be 
compared with results for a coarsening foam in Ref. 
where domains are defined physically by actual bubbles 
rather than by a Voronoi construction. For dynamics 
that are spatially uniform, the persistent area vanishes 
uniformly at the same rate throughout the whole sample, 
with only random uncorrelated variation between neigh- 
boring regions. This is not the case in Fig. |4j where 
the spatially-heterogeneous nature of the dynamics is ev- 
ident in the long swaths of white caused by chains of fast- 
moving beads. This is most pronounced for delay times 
between about 5 and 50 s, and corresponds closely with 
structure in the average velocity field seen for example in 

Fig. m 

The average persistent area and a dynamical suscepti- 
bility quantifying spatial heterogeneity may now be de- 
fined, in analogy with the previous section, by moments 
of A t (r) averaged over all times t: 



Air) = (A t (r)) 
X a{t) = N[{A t (rf)~{A t {r)f 



(6) 
(7) 



size. As an example, the decays of A t (r) vs r plotted in 
Fig. [5^i display considerable variation for different choices 
of starting time t. The resulting susceptibility, Xa{t), 
plotted underneath in Fig. [5)3, displays a peak as ex- 
pected near the characteristic decay time for the average 
persistent area. The height of this peak, x*ai an d the 
corresponding average persistent area A* , give the typi- 
cal number of beads in the dynamic heterogeneities as 



n A 



Xa 



(A± - A )(Ai - A*) ' 



(8) 



The variance is multiplied by the number N of beads 
in the sample, so that xa(t) is independent of system 



by the same arguments giving Eq. Here, Ao is the 
contribution to the persistent area by the beads in the 
fast-moving heterogeneities and A\ is the contribution 
from the remaining less mobile beads. Based on graphics 
like Fig. [4] we estimate Ao « and A\ ss 3/8. To remove 
this uncertainty, an alternative approach might be to en- 
force Aq — and A\ = 1 by defining an overlap order 
parameter such that the contribution from each bead is 
a step function that drops from one to zero when all its 
persistent area vanishes. 

The behavior of A(t) and xa{ t ) are displayed in Fig. [5] 
for different bead packing fractions, <j>, using the same se- 
quence of values and the same symbol/color codes as in 
Fig. [3] The decay time of A(t) and the peak position t* a 
of xa(t) coincide and become longer for increasing </> on 
approach to jamming. For all area fractions, the value of 
the persistent area is nearly constant, A* — 0.24 ± 0.04, 
when the susceptibility reaches its peak. The peak height 
X\ a l so increases as the dynamics slow down, indicative 
of a diverging dynamical correlation length. This behav- 
ior strongly parallels that in Fig. [3] based on an overlap 
order parameter with a particular step- function cutoff. In 
fact, the ratio of peak heights is nearly constant for all 
area fractions: (x*a/xD = 0.10 ± 0.02. The only striking 
qualitative difference is that A(r) exhibits a two-step de- 
cay at high (f>, whereas Q(r) does not. This is because the 
persistent area begins to decay at very short delay times, 
when the beads experience small-displacement ballistic 
motion, rattling within a "cage" of unchanging near- 
est neighbors. The shoulder in A(t) corresponds to the 
crossover from ballistic to subdiffusive motion observed 
in the mean-squared displacmcnt. For a typical cutoff I 
comparable to bead size, this so-called (3 relaxation is not 
detected by the traditional overlap order parameter Q(r), 
which decays only due to the cage-breakout a-relaxation 
process. Furthermore, no shoulder is evident in Q{1, r) 
for any choice of I in Fig. O Besides being naturally and 
uniquely defined by topology, the persistent area order 
parameter thus has the added advantage of being able to 
capture the two-step f3 and a relaxation processes char- 
acteristic of glass-forming systems. 



B. Persistent bond 

A second order parameter may be defined naturally 
from the set of nearest neighbors, or "bonds", specified 
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FIG. 4: Persistent area images vs delay time, r, as labeled, for air-fluidized beads at an area packing fraction of <f> = 0.792. The 
shrinking black regions represent persistent area, which has remained inside the same Voronoi cell since zero delay time. Each 
graphic represent the central 7.5 cm diameter region of the entire 17.7 cm diameter sample. The dynamics are most spatially 
heterogeneous at delay time t* a = 34 s where the susceptibility xa(t) reaches a maximum. 
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FIG. 5: (Color online) Time dependence of (a) the instanta- 
neous and average persistent areas, denned by Eqs. 1|5|6|) . re- 
spectively, and (b) the corresponding dynamic susceptibility, 
defined by Eq. ([7]), for air-fluidized beads at an area pack- 
ing fraction of cj> = 0.792. The starting times t for the in- 
stantaneous persistent areas are separated by 180 s, and thus 
represent statistically independent initial configurations. 



FIG. 6: (Color online) Time dependence of (a) the average 
persistent area, and (b) the corresponding dynamic suscep- 
tibility, for air-fluidized beads at different area fractions as 
labeled. 



by the Delaunay triangulation dual to the Voronoi con- 
struction. In particular, we define the instantaneous per- 
sistent bond number B t (r) as the fraction of all bonds at 
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FIG. 7: (Color online) Persistent bond images vs delay time, r, as labeled, for air-fluidized beads at an area packing fraction 
of cj> = 0.792. Each graphic represent the central 7.5 cm diameter region of the entire 17.7 cm diameter sample. The dynamics 
are most spatially heterogeneous at delay time t% = 180 s where the susceptibility xb (t) reaches a maximum 
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FIG. 8: (Color online) Time dependence of (a) the average 
persistent bond, and (b) the corresponding dynamic suscep- 
tibility, for air-fluidized beads at different area fractions as 
labeled. 



time t that remain unbroken across the interval t — *■ t+T. 
The progressive breaking of bonds throughout the sample 
is illustrated in Fig. [71 at the same area fraction as in the 
persistent area fraction diagrams of Fig. 2] By compari- 



son, the persistent bond number decays more slowly, due 
to neighboring beads that move a large distance while 
remaining next to one another. Also by comparison, the 
string-like swirls of the dynamic spatial heterogeneities 
are less evident by casual inspection of the persistent 
bond diagrams than for the persistent area diagrams. 
Nevertheless, bond lifetimes and the spatial heterogene- 
ity of broken bonds have been used to characterize sim- 
ulations of glassy systems [23[ . 

The average persistent bond number, the related dy- 
namic susceptibility, and the size of the dynamic hetero- 
geneities, respectively, are given in precise analogy with 
the previous sections as follows: 



B(r) 
Xb{t) 



N [(B t (T 
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(9) 
(10) 

(11) 



The trends in this order parameter and susceptibility are 
displayed vs r in Figs.[8^,-b for the same sequence of area 
fractions as in the analogous plots based on self-overlap 
and persistent area. Here, as before, the order parameter 
decays more slowly and the peak height increases as the 
area fraction is increased towards jamming. Similar to 
Q(t), but in contrast with A(t), the average persistent 
bond number B(t) exhibits a one-step decay. By inspec- 
tion, we estimate Bq ss 1/3 and B\ « 2/3. Note that 
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for all area fractions, the value of the persistent bond is 
nearly constant, B* = 0.46 ± 0.08, when the susceptibil- 
ity reaches its peak. Also, the ratio of peak susceptibility 
heights is nearly constant: (xh/xt) = 0-13 ± 0.05. 



IV. COMPARISONS 

Three dynamic susceptibilities, X4( r )i Xa(t), and 
Xb{t), have now been discussed for characterizing spa- 
tial dynamic heterogeneities in terms of the variance of 
order parameters based respectively on self-overlap, per- 
sistent area, and persistent bond. The nature of these 
order parameters, and their possible relation with local 
structure, are compared in Fig. [5] by snapshots of the 
system at eight different instants in times. Here the time 
increment is chosen as t a /2, half the time delay at which 
Xa(t) reaches its peak, so that there is a reasonable 
balance of continuity and evolution between successive 
snapshots. While each row represents a different time, 
the first three columns represent local structure and the 
last three represent local dynamics. In terms of struc- 
ture, Voronoi cells are shaded according to (a) their num- 
ber of sides; (b) their circularity shape factor, equal to 
perimeter-squared divided by An times area [2(| 0, H|| ; 
and (c) the reciprocal of their area, as a measure of the 
local density. In terms of dynamics, the average velocity 
vectors and the persistent areas are depicted at t* a and 
the persistent bonds are depicted at Tg, so that the spa- 
tial dynamic heterogeneities are maximally emphasized. 
By inspection, we note the following points. First, the 
spatial structure of the dynamics revealed by the aver- 
age velocity vectors and the persistent areas is clearly 
related. Both measures give a clear feeling for string-like 
swirls of neighboring beads that is intermittently excited 
against a background of less mobile beads. By contrast 
the persistent bonds give no such feeling of motion. And 
while there are more broken bonds near the fast mov- 
ing domain, the correspondence isn't one-to-one because 
neighboring mobile beads may remain nearest neighbors. 
Second, there is no apparent correlation between any of 
the measures of structure and the location of the spatial 
dynamic heterogeneities. It is not evident how to predict 
when or where a dynamic heterogeneity will occur based 
on usual measures of local structure. 

Next we compare the characteristic times and domain 
sizes for dynamic heterogeneities, as deduced from the 
three order parameters and susceptibilities. The growth 
of these scales as a function of increasing bead area frac- 
tion are displayed in Fig. [TUJ Note that the time scales 
are different for the three susceptibilities. As expected, 
Tg is the largest while t* a is only slightly larger than r| . 
However all three times appear to grow in proportion to 
one another, suggesting that they probe similar physics. 
Further reinforcing this point, the domain sizes are in- 
distinguishable for the three susceptibilities. Thus the 
overlap order parameter with discretionary choice of step- 
function cutoff at I — 0.5d s is not as arbitrary as it may 



seem, since it reveals the same picture of spatially het- 
erogeneous dynamics as the persistent areas and bonds, 
which arc uniquely defined by topology. 

Next we consider the functional form of the growth of 
dynamical heterogeneities on approach to jamming, ex- 
actly per Ref. |17|. Thus in Figs. [TUa-b. data for the 
peak times and domain sizes are shown on log-log plots 
vs |0 - ^d" 1 . With a single choice of <f) c = 0.79 ± 0.02, 
the data can all be well- fit to a power of \<j>— C | _1 where 
the exponent is v = 0.72 ± 0.04 for the domain sizes and 
is zv — 1.03±0.06 for the peak times. Such power-law di- 
vergences on approach to (f> c below random-close packing 
are consistent with a mode-coupling theory, where the 
mode-coupling temperature is above the glass-transition 
temperature. The value of the correlation length expo- 
nent v is consistent with prior simulation results based 
on finite-size scaling analysis of random close packing 
fraction (26[, on the size of the disturbance away from 
a perturbation [27] ] , and on scaling analysis of stress and 
strain rate [28|. Slightly smaller exponents for the dy- 
namical correlation length have been predicted for the 
approach 4> c from above [1^, [3(| [M| . The observed time 
scale exponent is less than that for simulation results 
based on stress relaxation , for the same model 

as in Ref. 28]. However the power-law description is 
not unique. As for the glass transition in molecular liq- 
uids, the growth of characteristic time and length scales 
can also be well-fit to Vogel-Tammann-Fulcher (VTF) 
equation, exp(_E/|0— (j) \). Thus in Figs. [TUb-d. data are 
shown on semi-log plots vs \<j) — 0o| _1 - With a single 
choice of 4>q — 0.83 ± 0.01, which corresponds to random 
close packing, the data follow the VFT equation where 
E = 0.20±0.01 for the peak times and E = 0.14±0.01 for 
the domain sizes. Altogether, as reported in Ref. [13], the 
growth of dynamical time and length scales on approach 
to point-J in the system of air-fluidized beads happens 
in close quantitative analogy to the glass transition in 
liquids. 

Before closing we make one final comparison of the 
four-point dynamic susceptibility for runs of different du- 
ration. In particular, we compute peak heights xX f° r 
subsets of duration T of the full 120 minute runs and 
normalize by the long-duration value xX- Results are 
plotted in Fig. [11] vs duration time, normalized by the 
long-duration peak time r\ , for all area fractions consid- 
ered above. Note that the normalization causes all data 
to collapse, even though peak heights and times grow on 
approach to jamming. Also note that the long-duration 
value of the peak height is attained only if the run is 
about ten times longer than the peak time r| . Intuitively, 
if a run is not very long compared to the decay time of 
the overlap order parameter, then too few configurations 
are sampled and the variance is underestimated. This 
effect will limit the proximity to which point-J may be 
approached. For example, if run duration is held fixed 
and the area fraction is gradually increased, then the 
peak height and the corresponding dynamical correla- 
tion length will initially grow but will eventually appear 
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circularity 
x = 



local density ave. vel. vectors persistent area persistent bond 

T — T^* T — T^* T — T Tg* 




FIG. 9: (Color online) Graphic illustration of three measures of structure and three measures of dynamics, at eight instants of 
time. Column- 1: Voronoi cells shaded according to their number of sides Z; cells with many sides are darker than cells with 
few. Column-2: Voronoi cells shaded according to their circularity f , defined by perimeter squared divided by 4n times area; 
more circular cells are darker than less circular. Column-3: Voronoi cells shaded according to local density, proportional to the 
reciprocal of their area and averaged over the time interval r%\ cells with low local density are colored darker than cells with 
high local density. The standard deviation of density values equals about one-tenth the average density. Columns 4-6: average 
velocity vector fields, persistent areas, and persistent bonds, respectively, at delay times as labeled. 
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FIG. 10: (Color online) Growth with increasing area fraction of (a,c) times at which the three dynamic susceptibilities Xq( t )> 
Xa{t), and xb(j) become maximum and (b,d) the domain sizes for the spatial dynamic heterogeneities as deduced from Eqs. (j4] 
HlQjJ. Best fits to power law divergence as \<f> — cj> c \ yield a critical density of 4> c = 0.79 ± 0.01 and exponents as labeled. 
Best fits to exponential divergence as exp(_E/|0 — 4>o\) yield a Vogel-Thammann-Fulcher critical density at <f>o ~ 0.83 ± 0.01 and 
energies as labeled. Solid symbols are for the 120 minute runs featured in previous plots; open symbols are for 20 minute runs. 



to decrease as r| grows beyond about one-tenth of the 
run duration. This can be seen in the inset of Fig. [TTJ To 
eliminate such erroneous finite-time artifacts, the experi- 
mental observation time must be at least a decade longer 
than the time it takes for the overlap order parameter to 
vanish. Conceivably, an alternative might be to take a 
smaller cutoff in the self-overlap order parameter or to 
alter the proportionality constant relating susceptibility 
and domain size according to the experimental condi- 
tions. 



V. CONCLUSION 

Spatially-heterogeneous dynamics may be character- 
ized successfully using four-point susceptibilities based 
on many choices for the configurational order parame- 
ter. For a self-overlap order parameter with the seem- 
ingly ad-hoc choice of a step- function cutoff at d s /2, the 
characteristic dynamical correlation length for the size 
of intermittent mobile domains is indistinguishable from 
those based on persistent areas and bonds in a unique 
topological description. This agreement requires a phys- 
ical picture for how to deduce domain size from the peak 
of the susceptibility, as in Eqs. (|4I8I11|) . Nevertheless, 
some differences still exist in the detailed form of the 
three susceptibilities. For example the peak times all dif- 
fer, though they always remain in constant proportion to 
one another. At early times, below the peak, the persis- 
tent area is the first to begin its decay since self-overlap 
and persistent bonds are insensitive to the small ballis- 
tic motion of grains as they rattle in a cages of fixed 
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FIG. 11: (Color online) Peak height xt vs observation time T, 
for different area fractions as labeled. Heights and times are 
normalized, respectively, by the actual long-duration values 
Xt and the peak times t%. To obtain good configurational 
averaging and accurate experimental results, the duration of 
the observation time must be about a decade longer than the 
decay time of the overlap order parameter. Inset: Apparent 
peak height vs area fraction for short and long observation 
times, as labeled. Symbols in the main plot denote different 
area fractions, as specified by corresponding symbols in the 
inset. 



nearest neighbors. At late times, past the peak, the per- 
sistent bond is the last to decay since a few neighboring 
beads remain in close contact as they diffuse through the 
sample. Thus, while all three susceptibilities give the 
same picture of the spatially heterogeneous dynamics, 
they offer different possibilities for characterizing short 
and long time dynamics. As applied to the change in 
behavior with increasing bead packing fraction, all three 
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give characteristic time and length scales for the dynam- 
ical heterogeneities that appear to diverge in accord with 
simulation results for supercooled liquids and for dense 
athermal systems of soft repulsive particles. 
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